Finite element approach for simulating quantum electron dynamics in a magnetic field 
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A fast and stable numerical method is formulated to compute the time evolution of a wave function 
in a magnetic field by solving the time-dependent Schrodinger equation. This computational method 
is based on the finite element method in real space to improved accuracy without any increase of 
computational cost. This method is also based on Suzuki's exponential product theory to afford an 
efficient way to manage the TD-Schrodinger equation with a vector potential. Applying this method 
to some simple electron dynamics, we have confirmed its efficiency and accuracy. 
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I. INTRODUCTION 



II. FORMULATION 
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Conventionally, wave functions have been represented 
as a linear combination of plane waves or atomic orbitals 
in the calculations of the electronic states or their time 
evolution. However, these representations entail high 
computational cost to calculate the matrix elements for 
these bases. The plane wave bases set is not suitable for 
localized orbitals, and the atomic orbital bases set is not 
suitable for spreading waves. 

To overcome those problems, some numerical methods 
adopted real-space representation to solve the time de- 
pendent Schrodinger equation In those methods, 
a wavefunction is descritized by grid points in real space 
and the spatial differential operator is approximated by 
the finite difference method (FDM). With those meth- 
ods, some dynamic electron phenomena were simulated 
successfully [0-^]. 

In the previous work |TT| ], we have formulated a new 
computational method for the TD-Schrodinger equation 
by using some computational techniques such as, the 
FDM, Suzuki's exponential product theory fl^-|l7|] , Cay- 
ley's form H and Adhesive operator. This method af- 
forded high-stability and low computational cost. 

In the field of engineering, for example, numerical anal- 
ysis of fluid dynamics or of strength of macroscopic con- 
structions, the finite element method (FEM) has been 
widely and traditionally used for approximating the ap- 
propriate partial differential equations. Recently, the 
FEM has been found useful for the time-independent 
Schrodinger equation of electrons in solid or liquid mate- 
rials (To). 

In this paper, we have utilized the FEM for solving 
the TD-Schrodinger equation as an extension of the pre- 
vious work JTl[ . By using Cayley's form and the FEM, 
this method affords high-accuracy without any increase 
of computational cost. Moreover, we have formulated a 
new efficient method which manages the time evolution 
of a wave function in a vector potential or in a magnetic 
field. These techniques are especially useful for simu- 
lating dynamics of electrons in a variety of meso-scopic 
systems. 



In this section, we formulate a new method derived by 
the FEM and a new scheme to manage a vector poten- 
tial efficiently. Throughout this paper, we often use the 
atomic unit h — 1, m = 1, e = 1. 



A. FEM for the TD-Schrodinger equation 

First, we utilize the FEM for the time evolution of 
a wave function in a one-dimensional closed system de- 
scribed by the following TD-Schrodinger equation: 



dt 



2m dx 2 



(1) 



The FEM starts by smoothing the wavefunction 
around a grid point. We smoothed ijj(x) around a grid 
point Xi by eq. (||), as illustrated in Fig. [l| 

^{x,t)=^ i {t)+^ l i {t){x-x i ) + -^l[t){x-x i f , (2) 



where 



ipi(t) = 1p(Xi,t) , 

i>i+i{t) - ipi-i(t) 



2Ax 

&+i(t)-2^(t)+V<-i(t) 



(3) 




FIG. 1. The FEM starts by smoothing the wavefunction 
around a grid point. The wavefunction is supplemented by a 
quadratic equation. 
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By substituting eq. (||) for eq. (0), tp(x, t) is expressed Using these notations, eq. (0) is expressed simply as 



ip(x,t) = u a ,i(x)ip i ^i(t) + u oi (x) ipi(t) +uu(x) 1pi+l(t) 



(4) 



where u a i(x), u Q i(x) and u\,i(x) are the base functions 
defined below: 



u ai (x) 



u m (x) = 1 - 
u hi (x) 



(x Xj) (x x^ 
2Ax 2 2Ax 

(x-x,i) 2 



Ax 2 ' 

(x-x t ) 2 {X-Xi) 



(•5) 



2Ax 2 2Ax 



Substituting eq. (Q) for eq. ([!]) and multiplying both 
side of the equation by the base function u i{x) and in- 
tegrating by x in the range [a^_i, Xi + {\ as 



ih / dxu oi (x) u oi (x) ipi(t) 

J Xi—l 

+ u ai (x) ipi+i (*) + u hi (x) ipi- 1 (t) 

f 1 r 

da;woi(a;) dlu i{x) ipi(t) 



hf_ 

2m 
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+ diu ai (x)tp i+1 (t) + d^u bi (x)if)i-i(t) , (6) 



the following formula is obtained after some algebra: 
ih-^ [ipi-i{t) + 8^i(t) + ^<+i(t)] 



2mAx 2 



[^-i(t)-2^i(t)+Vi+i(*)] • (7) 



To simplify the expression, it is useful to define a vector 
and two matrices as below: 



(8) 



8 10 

1 '■. '•. 

'■. '■• 1 

18 



D 



-2100 
1 '•• '•• 

'■. '■• 1 

1-2 



(9) 



Clearly S and D satisfy the following equation: 

1 



S = I 



10 



-D 



(10) 



iTiS 



(11) 



dt 2mAx 2 

Equation (^) is the finite element equation for this case. 

It has been thought that the existence of the matrix 
S is troublesome since the inverse of this matrix is re- 
quired to obtain the time derivative of the wave function, 
namely, 



(12) 



dt 2mAx 2 
However, we have found that this differential equation 
is easily solved by using an approximation called Cayley's 
form. The formal solution of eq. (|l2| ) is given by, 

ih At 



ip(t + At) = exp 



-S *D 



m ■ (13) 



-2m Ax 2 

The exponential operator is approximated by Cayley's 
form: 



i/j(t + At) 



ih At 
Am Ax 2 



-S D 



I 



ih At 



(14) 



-S D 



Am Ax 2 

Multiplying both the numerator and the denominator of 
the righthand side by the matrix S and using the relation 
(0), the required formula is obtained: 

At 



tp(t + At) 



ih 



Am c tf Ax 2 



D 



ih At 

Am* s Ax 2 



(15) 



D 



where m c s is an "effective mass" of an electron defined 
as 



.2 Ax 2 



l- 



5 At 



(16) 



In this way, the solution of the partial differential equa- 
tion, eq. (|I]) is computed by eq. (Il|) with the concept of 
the FEM. It is quite a remarkable result that formula 
eq. (1l5[) is almost the same as the formula derived by the 
FDM ||ll|. In this time evolution, the norm of the wave 
function is exactly conserved since the time evolution op- 
erator appearing in eq. (|l|) is strictly unitary. Moreover, 
accuracy is dramatically improved without any increase 
in the computational cost, as demonstrated in the next 
section. 

It is easy to extend this idea for two-dimensional 
systems, since the time evolution operator in a two- 
dimensional system is decomposed into a product of the 
time evolution operators in one-dimensional systems [ [il"| . 
The approximated solution utilizing the FEM is given by 

ih At 



ip(r, t + At) 



4m c ff Ax 2 



D, 



ih At n 
D„ 

4m Gff Ay 2 



ih At 

Am* s Ax 2 



D, 



ih At 



Am* Ay 



(17) 
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where and T> y are the finite difference matrices along 
the x and y axes respectively, and their appearances are 
the same as D defined in eq. (H). 



Moreover, we have found that the hybrid decomposi- 
tion JItJ is rather easy in this case. Note the following 
identity: 



B. Evolution in a magnetic field 

Though there are many interesting phenomena in a 
magnetic field, there has been no efficient methods that 
numerically manage the dynamics in a magnetic field as 
far as we know. We have improved our method to afford 
an efficient way to solve the TD-Schrodinger equation 
with a vector potential given as below 



at 



2m 



ti 



A ^(r,f) 



(18) 



In this subsection, we present the method for only the 
case of a two-dimensional system lying on the xy plane 
subjected to a uniform external magnetic field along the z 
axis. We do not mention the case of a non-uniform mag- 
netic field specifically, but the extension of the method is 
straightforward. We adopt the following vector potential 
A for this magnetic field: 

A = {-By, 0, 0) T . (19) 

The TD-Schrodinger equation of this system is given by 

dtp(r,t) 



dt 



(20) 



The strict, analytical solution is also given by an expo- 
nential operator: 



V»(r,t + At) = exp\^-At(d x -^ByY + ^Atd, 
12m V n > 2m 



[dl[(d x -iayr,d 2 y ]]^~8a^ 



(24) 



Then, equation ( pl|) is approximated by the following 
fourth-order hybrid exponential product: 



tp(r,t+ At) = exp 



At—[- 

2m V6 



le 



x exp y-\-—Bxyj exp 



x exp 



x exp 



x exp 



2m V6 



Bxy ) exp 
I e 2 B 2 At 2 



[At ifi 

~~2 2m 
[2 At ih 2 

~2^i v 

At in 2 



e 2 B 2 At 2 
72m 2 c 2 

exp \ ——Bxy 



2 2m d 



exp 



(->») 



72m 2 c 2 



V-(r,i)+0(At 5 ; 



(25) 



The exponential of the magnetic field just changes the 
phase of the wave function, so it is very easy to compute. 
Therefore, this method is adaptable to systems subjected 
to a magnetic field. The outline of the procedure for a 
two-dimensional system subjected to a magnetic field is 
schematically described by Fig. g. 
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FIG. 2. The procedure for a two-dimensional system sub- 
jected to a magnetic field. Here Bxy means the operation of 
. . the exponential of the magnetic field. In this way, the phase 
' of the wavefunction is turned forward before the operation of 
(21) Cayley's form along the x-axis and is turned backward after 
Cayley's form. 



Note the following identity: 



exp 



Af 



m 

'2m 



(d x 



ie n 2 

n By 



x exp 



At^d 2 
2m 



— exp 



exp 



(22) 



III. APPLICATIONS 



A. Comparison between FDM and FEM 



Equation ( f2l| ) is approximated by the following second- 
order exponential product: 



ip(r,t + At) = exp 



At in 2 
Y2^ v 



exp i 



x exp 



2m 



x exp 



At in 2 
Y2^ y 



(+-B X y 
exp (--^Bxy) 
i>(r,t) + 0(At 3 ) 



(23) 



In this subsection, we briefly compare Cayley's form 
and other conventional methods by simply simulating a 
Gaussian wave packet moving in a one-dimensional free 
system as illustrated in Fig. 0. 
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FIG. 3. The model system for comparison with the con- 
ventional methods. 256 computational grid points are allo- 
cated in the physical length 8.0a.u. A Gaussian wave packet 
is placed in the system, whose initial average location x a and 
momentum p are set as x = 2.0a.u. and p — 12.0a.u., re- 
spectively. 

The TD-Schrodinger equation of this system is simply 
given by 



,dtp(x,t) 
1 dt 



d 2 

-f V>(M). 



(26) 



The wavefunction at the initial state is set as a Gaussian: 



ip(x,t = 0) = 



1 



■. exp 



AW 2 



ip x 



, (27) 



where W — 0.25a.u. x — 2.0a.u. p — 12.0a.u. The evo- 
lution of this Gaussian is analytically derived as 



ip(x,t) 



x exp 



1 



^/2nW 2 + (ir/2)(t/W) 2 

(x — X — p t) 2 



4w 2 + (t/wy 



ip X 



(28) 



Therefore, the average location of the Gaussian (x(t)) is 
derived as if it is a classical particle: 



(x(t)) = (x(t = 0))+p o t 



(29) 



This characteristic is useful to check the accuracy of the 
simulation. 

Cayley's form with the FDM is given by 

where d 2 is approximated by a finite difference matrix as 



(31) 



Meanwhile, Cayley's form with the FEM is given by 

(32) 
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"left + iAt/Ad 2 , . , 



& - iAt/AO 2 

where the spatial differential operator is approximated in 
the ordinary way and m e s is the effective mass: 



f 



1 2 Ax 2 



5 At 



(33) 



d 2 is approximated by eq. (plh . 



We have simulated the motion of the Gaussian by those 
methods. Figure ^ shows the error in the average momen- 
tum. The errors are evaluated in the following way: 



e(At/Ax*)= {x{t = Z )) - X ° - 1) „ 



T 



(34) 



N-l 



(x(t)) = Ax ^ x i\Mt)\ 2 in FDM. (35) 



Ax 



i=0 
N-l 



x(t)) = —Re Xi^*(24tpi + 4^+i + 4^_i 



-4>i+2 - ipi-2) in FEM. 



(36) 
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FIG. 4. The errors in the average momentum computed 
by Cayley's form with the FDM and Cayley's form with the 
FEM. The error of the FEM is smaller than that of the FDM. 
The spatial slice is set as Ax — l/32a.u. 

It is found that the accuracy is dramatically improved 
by using the FEM. It is remarkable that in spite of the 
improvement of accuracy, the computational cost does 
not increase at all. 



B. Cyclotron motion 

We demonstrate the cyclotron motion in the frame- 
work of quantum mechanics. We have simulated the mo- 
tion of a Gaussian wave packet in a uniform magnetic 
force as illustrated in Fig. ||. 
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FIG. 5. The model system for the cyclotron motion. This 
system is subjected to a static magnetic force perpendicularly, 
and it is surrounded by infinitely high potentials. 64 x 64 
computational grid points are allocated in the physical length 
8a.u. x 8a. u. The strength of the static magnetic force B is 
set as 2a.u. A Gaussian is placed as the initial state of the 
wavefunction, whose average location and momentum are set 
as (6a.u., 4a.u.) and (Oa.u., 4a.u.), respectively. The time slice 
is set as At = l/64a.u. 



The initial wavefunction ip(r,t = 0) is set as the fol- 
lowing Gaussian: 




t=o 



t=3/8 



t=6/8 



t=9/8 



V>(r, t = 0) 



1 



V2irW 2 





r |r-r | 2 " 




exp 


AW 2 - 


exp 



ieB 



(x-L/2)y 
(37) 



where r Q is set as x Q = 6a. u., y = 4a. u. and W is set as 
0.5a.u. 

The initial density p(r, t = 0) and the initial current 
density j(r, t — 0) derived from this wave function are as 
follows: 




t=12/8 



t=15/8 



t=18/8 



t=21/8 



p(r,t = 0) 
j(r,* = 0) = 



2irW- 



exp 



r - r D 



2W 2 



2 Bp(r) 



mc 



We adopt a gauge of the vector potential A as 



A = (-By, 0, 0) 1 



(38) 



(0, x- L/2, 0) T . (39) 



(40) 



FIG. 6. The evolution of the density and the current vec- 
tor. The Gaussian is observed to circle around. 
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In classical mechanics, the average momentum of this 
Gaussian at the initial state is evaluated as 



Po = -\(S)\ = —\x -L/2\ 



(41) 



FIG. 7. The orbit of the average location of the wave 
packet. The radius of this circular trace is estimated as 2. Oa.u. 
The initial average location and momentum of this Gaussian 
are set as (6a.u., 4a.u.) and (Oa.u., 4a.u.), respectively. 



This means the classical cyclotron radius is \x — L/2\. 

Some snapshots of the simulation time span are illus- 
trated in Fig. § The average location of the wave packet 
is observed to circle around as plotted in Fig. 0. 



This trace is not a perfect circle but a swirl due to the 
reflection by the closed walls around the system. 

A more perfect circular trace is observed by enlarging 
the system or shortening the cyclotron radius to reduce 
the effect of the reflection. Figure || shows the result of 
another simulation. 
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FIG. 8. The another orbit of the average location of the 
wave packet. The radius of this circular trace is estimated 
as l.Oa.u. The initial average location and momentum of this 
Gaussian are set as (5a.u., 4a.u.) and (Oa.u., 2a.u.), respec- 
tively. 

These results afford good agreement with the result by 
classical mechanics. 



where d and D mean the width of the slits and the span 
of the slits respectively. Thus D — d is the length of the 
wall where a magnetic flux goes through. 

In an analogy to semi-classical photon interference, the 
electron interference pattern I(x) in this AB system is 
approximately described by the following form: 



I(x) oc 



21 . 


\kd ~ 




\kD 








cos 


hi*- 




kdx 


V2l X . 




2h- 



(44) 



In the above, £ is a y coordinate where the pattern is 
evaluated. 

Figure [if] shows the result of this simulation for the 
case of no magnetic flux, $ = 0. These data were taken 
soon after the pattern appeared in order to prevent the 
pattern from extra interference due to the reflected waves 
from side walls. The interference pattern basically agrees 
with the semi-classical one derived from eq. (^J) . 



C. Aharonov-Bohm effect 



We demonstrate Aharonov-Bohm effect by simulating 
an electron dynamics on a system as illustrated in Fig. |^. 
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FIG. 9. The model system for the Aharonov-Bohm effect. 
The shape of this system is rectangular. A double-slit lies 
at the center. A magnetic flux <J? goes through a wall lying 
between the slits. 64 x 128 computational grid points are 
allocated in the physical size 8a. u. x 16a.u. The initial wave- 
function is set as a plane wave k in front of the double-slit. 
The time slice is set as At = I/64a.u. 

The vector potential is constructed as follows: 
A(x,y) = (0, Ay(x), Of ; A v (x) = - f dx'B(x',y) . 

J-L/2 



(42) 



Thus A y (x) has a finite value only inside the right slit 



—B(D — d) : inside the upper slit. , .„.. 
: in other area. ' 









FIG. 10. The interference pattern observed in the back of 
the double-slit and at the line y = £ = L/4 in a case of no 
magnetic flux, $ = 0. The solid line indicates the numerical 
result; the dashed line indicates the semi-classical one derived 
from eq. Q) . 

Further, the results for the case of magnetic flux 
$ = h/2e and $ = h/e are shown in Figs. |ll| and [l^, 
respectively. The patterns are observed to shift to the 
right-hand side, and these behaviors also agree with the 
semi-classical one. However, the patterns are different 
from the the semi-classical one in their details. This is of 
course due to the quantum effect. 
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FIG. 11. The interference pattern observed in the back 
of the double-slit and at the line y = I = L/4 in a case 
of <E> = h/2e. The solid line indicates the numerical result; 
the dashed line indicates the semi-classical one derived from 
eq. @. 
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FIG. 12. The interference pattern observed in the back 
of the double-slit and at the line y = I = L/4 in a case 
of $ = h/e. The solid line indicates the numerical result; 
the dashed line indicates the semi-classical one derived from 
eq. @. 



IV. CONCLUSION 

We have improved the computational method for the 
time-dependent Schrodinger equation by utilizing the fi- 
nite element method and by formulating a new scheme 
for a magnetic field. We have found that by using the 
FEM, the accuracy of the simulation is dramatically im- 
proved without any increase in the computational cost. 
We have also found that the new scheme is quite efficient 
for simulating systems in a magnetic field. 

This computational method is especially useful for sim- 
ulating dynamics of electrons in a variety of meso-scopic 
structures. 
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